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Abstract.  Flow  of  two  immiscible,  incompressible  fluids  in  a  porous  medium  is  typ¬ 
ically  described  by  a  nonlinear  advection-diffusion  equation  for  one  of  the  fluid  satura¬ 
tions.  The  diffusion  coefficient,  which  represents  the  effect  of  capillary  forces  on  the  flu¬ 
ids,  is  zero  when  the  medium  is  locally  saturated  by  either  fluid  since  in  these  limiting 
cases  the  effects  of  capillary  forces  tend  to  zero.  This  degeneracy  in  the  second-order  term 
usually  gives  rise  to  the  qualitative  property  that  perturbations  in  saturation  propagate 
with  finite  speed  through  regions  that  are  fully  saturated  by  either  fluid.  This  qualita¬ 
tive  property  is  physically  realistic.  In  this  work  we  show  that,  under  certain  choices  of 
constitutive  relations  and  modeling  approximations,  the  finite  speed  of  propagation  prop¬ 
erty  is  lost,  despite  the  fact  that  the  diffusion  coefficient  is  degenerate.  The  loss  of  fi¬ 
nite  speed  of  propagation  is  due  to  unbounded  derivatives  in  the  closure  relations  as  the 
medium  becomes  saturated  by  wetting  phase.  We  present  analytical  and  numerical  so¬ 
lutions,  compare  solution  dynamics  that  display  finite  and  infinite  speed  of  propagation, 
and  provide  a  brief  account  of  numerical  difficulties  related  to  the  degenerate  coefficients. 


1.  Introduction 

Models  of  two-phase  flow  in  porous  media  are  fundamen¬ 
tal  tools  for  understanding  subsurface  water  resources.  They 
describe  soil  moisture  in  the  unsaturated  zone  as  well  as 
transport  of  industrial  solvents  in  aquifers.  Models  of  two- 
phase  flow  must  mimic  static  and  dynamic  conditions  in 
two-phase  systems  to  be  useful,  whether  they  are  used  to 
make  predictions,  to  design  remediation  strategies,  or  simply 
to  provide  physical  intuition.  Nevertheless,  in  representing 
nature  as  a  mathematical  model,  some  natural  behavior  is 
lost  and  some  unnatural  behavior  is  gained.  In  this  respect, 
good  models  present  acceptable  trade-offs  between  physi¬ 
cal  and  non-physical  behavior.  This  technical  note  focuses 
on  the  way  perturbations  propagate  in  solutions  of  certain 
two-phase  model  equations  and  the  consequences  of  this  be¬ 
havior  on  numerical  methods.  The  more  important  question 
of  whether  the  solution  behavior  is  an  acceptable  approxi¬ 
mation  of  nature  would  need  to  be  addressed  by  experiments 
or  deeper  theoretical  results. 

We  consider  models  based  on  continuum  mass  balance 
equations  and  the  multi-phase  extension  of  Darcy’s  law,  de¬ 
rived,  for  example,  using  the  representative  elementary  vol¬ 
ume  approach  [Bear,  1972].  The  fluids  are  assumed  to  be  in¬ 
compressible,  and  the  media  is  assumed  to  be  incompressible 
and  immobile.  This  approach  results  in  two  conservation 
laws  relating  the  saturations  of  the  fluids  to  the  fluid  pres¬ 
sure  gradients  and  to  the  gravitational  forces  via  nonlinear 
proportionality  constants  called  the  relative  permeabilities. 
Mathematical  closure  is  obtained  by  following  the  approach 
outlined  in  Parker  et  al.  [1987],  which  uses  the  capillary 
pressure  model  described  in  van  Genuchten  [1980]  and  the 
permeability  model  described  in  Mualem  [1976].  The  result 
is  a  coupled  system  of  second-order  quasilinear  equations 
that  are  of  degenerate  parabolic  type.  That  is,  except  at 
two  discrete  values  of  saturation,  the  diffusion  coefficients 
are  strictly  positive,  and  all  of  the  nonlinear  coefficients  are 
at  least  once  differentiable.  In  one  dimension  the  system 
of  equations  can  be  reduced  further  to  a  single  degenerate 
parabolic  equation.  This  work  focuses  strictly  on  the  behav¬ 
ior  of  the  one-dimensional  model  equation. 


The  qualitative  behavior  of  solutions  of  partial  differential 
equations  is  particularly  important  in  the  case  of  spatially 
variable  or  nonlinear  coefficients,  where  analytical  solutions 
are  not  generally  available.  When  something  is  known  about 
the  qualitative  behavior  of  solutions,  such  as  a  maximum 
principle  or  conservation  of  mass,  these  properties  can  be  re¬ 
covered  in  a  discrete  sense  by  carefully  designing  numerical 
methods.  Gilding  and  Kersner  [1996]  gives  a  very  complete 
classification  of  scalar  nonlinear  advection-diffusion-reaction 
equations  in  terms  of  finite  speed  of  propagation,  which  is  a 
qualitative  property  of  solutions. 

To  define  the  finite  speed  of  propagation  property  we  be¬ 
gin  by  considering  non-negative  functions  u(x,  t)  (i.e.  solu¬ 
tions  of  partial  differential  equations)  which  are  defined  on 
a  half  strip  (0,  oo)  x  (0,  T],  We  define  the  support  of  u  as 

P[t]  =  {x  £  (0,  oo)|u(x,  t)  >  0}  (1) 

and  the  front  or  free  boundary  as 

C  (t)  =  sup  P[t ]  (2) 

xG(0,oo) 

We  say  that  u  has  the  finite  speed  of  propagation  property 
if 

C(t)  <00  t  €  [0,  t)  (3) 

for  some  r  >  0.  That  is,  if  the  support  of  u  is  bounded  ini¬ 
tially  then  it  stays  bounded  over  finite  time  intervals.  The 
property  is  related  to  the  existence  of  finite  or  bounded  trav¬ 
eling  wave  solutions  of  the  partial  differential  equations,  as 
is  shown  in  Gilding  and  Kersner  [1996].  If  u  is  some  mea¬ 
sure  of  the  amount  of  a  substance  at  location  x  at  time  t, 
and  that  substance  is  merely  transported  by  natural  pro¬ 
cesses,  then  it  seems  reasonable  to  require  that  solutions  of 
the  model  equation  not  cause  the  substance  to  occupy  an 
unbounded  region  in  finite  time.  Equivalently,  we  expect 
that  £(t)  is  always  finite.  If  a  column  of  porous  media  is 
initially  saturated  with  fluid  W  and  another  fluid  N  is  in¬ 
jected  at  the  bottom  of  the  column,  then  the  finite  speed  of 
propagation  implies  that  N  will  not  flow  out  of  the  top  of 
the  column  instantaneously. 
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The  linear  advection-diffusion  equation  does  not  possess 
the  finite  speed  of  propagation  property:  solutions  with 
bounded  initial  support  are  instantaneously  positive  on  the 
entire  half  line.  This  behavior  in  the  mathematical  model  is 
acceptable  in  most  cases  because  the  solution  decays  rapidly 
enough  that  it  provides  a  good  approximation  of  nature. 
Nonlinear  hyperbolic  partial  equations  usually  have  finite 
speed  of  propagation.  On  the  other  hand,  nonlinear  hy¬ 
perbolic  equations  have  discontinuous  solutions,  which,  al¬ 
though  non-physical  in  many  applications,  are  often  good 
approximations  at  the  scale  of  observation.  The  objective 
of  this  work  is  to  demonstrate  a  non-physical  aspect  of  so¬ 
lutions  of  the  model  equation,  which  may  or  may  not  be  an 
acceptable  approximation  to  nature,  and  to  examine  how 
this  aspect  affects  numerical  methods. 

In  section  2  we  present  the  two-phase  flow  model  and 
Richards’  equation  for  air/water  systems.  We  also  derive 
some  asymptotic  approximations  to  the  nonlinear  parame¬ 
ters  that  will  be  used  to  formulate  simplified  model  equa¬ 
tions  and  to  analyze  the  speed  of  propagation.  We  apply 
the  results  of  Gilding  and  Kersner  [1996]  to  the  two-phase 
models  in  section  4  and  identify  the  range  of  parameters 
that  lead  to  infinite  speed  of  propagation.  In  section  4  we 
present  analytical  solutions  of  the  equations  for  several  spe¬ 
cial  cases  in  order  to  contrast  the  behavior  of  solutions  with 
and  without  finite  speed  of  propagation.  Finally,  in  section 
5,  we  describe  some  simple  numerical  approaches  for  solving 
the  equations  and  demonstrate  additional  consequences  of 
the  degenerate  form  of  the  equations. 

2.  Two-Phase  Flow  Models 

In  one  space  dimension,  the  standard  model  for  incom¬ 
pressible  two-phase  flow  in  a  homogeneous  porous  medium 
reduces  to  an  equation  for  the  wetting  phase  saturation  s 
(c.f.  [ Chavent  and  Jaffre,  1986]): 

St  +  [f(s)  -  d(s)sx}x  =  0  (4) 


where 


d 

f 


krw  krf 


dpc 


(-0  flw  ^  krw  +  krn  ds 


ki  fin 
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(5) 

(6) 


The  resulting  equation  is  the  saturation  form  of  Richards’ 
equation.  Since  the  density  and  viscosity  of  air  are  in  fact 
much  smaller  than  those  of  water,  the  assumption  that 
pn/ p-w  =  0  and  pn/pw  =  0  seems  like  a  reasonable  approxi¬ 
mation. 

2.2.  Effective  Saturation 

Most  closure  relations  for  two-phase  flow  are  based  on  a 
reduced  or  effective  saturation  defined  by 


where  sr  is  the  residual  saturation  and  ss  is  the  maximum 
wetting  phase  saturation.  This  reflects  the  fact  that,  if  only 
flow  processes  are  involved,  it  is  usually  impossible  for  either 
of  the  two  fluids  to  be  completely  removed  from  the  medium 
because  of  entrapment  mechanisms.  For  this  work  we  are  in¬ 
terested  in  what  happens  when  s  =  ss-  For  simplicity,  we 
take  s3  —  1  and  sr  =  0  so  that  se  =  s. 

2.3.  Relative  Permeability 

We  will  use  the  wetting  phase  relative  permeability  model 
derived  in  Mualem  [1976],  which  was  later  extended  to  the 
non- wetting  phase  in  Parker  et  al.  [1987]: 


The  analytical  form  of  these  functions  depends  on  the  cap¬ 
illary  pressure. 

2.4.  Capillary  Pressure 

The  capillary  pressure  is  expressed  as  a  function  of  satu¬ 
ration  according  to  the  model  presented  in  van  Genuchten 
[1980]: 


and  where  ki  is  intrinsic  permeability,  u  is  porosity,  pw  is 
wetting  phase  viscosity,  pn  is  non-wetting  phase  viscosity, 
pw  is  wetting  phase  density,  pn  is  non-wetting  phase  density, 
qt  is  total  fluid  velocity,  and  gx  is  gravitational  acceleration. 
The  nonlinear  constitutive  relations  are  the  relative  perme¬ 
abilities,  krw  and  krn,  and  the  capillary  pressure  pc,  which 
we  will  specify  after  considering  a  related  two-phase  model. 
We  will  refer  to  d  as  the  diffusion  coefficient  and  /  as  the 
advective  flux. 

2.1.  Richards’  Equation  for  Air /Water  Systems 


1  /  — l/m  ,  \  /  1 

Pc  =  ~  (se  '  ~  1)  (12) 

O'  ’ 

where  a,  n,  and  m  are  constants. 

2.5.  Van  Genuchten-Mualem  Closure  Relations 

Following  [ van  Genuchten ,  1980],  we  require  that  m  = 
1  —  1/n,  so  that  the  integral  in  the  permeability  relations 
above  can  be  integrated  exactly.  The  resulting  closure  rela¬ 
tions  can  be  written  as 


Our  intention  in  this  section  is  not  to  derive  Richards’ 
equation  rigorously  but  rather  to  clarify  how  commonly  used 
forms  of  Richards’  equation  relate  to  the  two-phase  model 
above.  The  original  derivation  of  Richards’  equation  was 
based  on  physical  principles  rather  than  as  an  approxima¬ 
tion  to  the  two-phase  flow  model  [ Richards ,  1931].  Starting 
from  the  two-phase  flow  equation,  if  we  assume  pn/ pw  =  0 
and  p-n/p-w  =  0  we  obtain 


/  = 


ki  kr 


d  =  - 


Pw  Qx 

LO  flw 

kikrw  dpc 


(7) 

(8) 


dpc 

dse 

krw 

krn 


=  si/2[i~(i-si/mr]2 

=  (l-Se)1/2(l-Se/m)2m 


(13) 

(14) 

(15) 


for  se  £  [0, 1]  where  0  <  m  <  1  and  a  >  0  are  constants. 

For  notational  simplicity  we  will  set  a  =  u>  =  ki  =  pw  = 
Pw  =  Pn  =  1,  and  pn  =  0.  Additionally,  we  consider  only 
cases  with  qt  =  0  and  gx  =  —  1  or  0.  The  case  with  gx  =  —  1 
is  called  gravity  segregation  or  counter-current  flow,  since 


copw  ds 
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Figure  1.  Advective  flux,  /,  versus  saturation,  s,  for 
two-phase  flow.  For  all  m,  /( 0)  =  /( 1)  =  0.  For  small 
m,  4^(1)  =  oo .  Consequently,  f  is  not  Lipschitz  contin¬ 
uous  at  s  =  1  for  small  m. 


S 


Figure  2.  Diffusion  coefficient,  d,  versus  saturation,  s, 
for  two-phase  flow.  For  all  m,  d(0)  =  d(l)  =  0.  For 
small  m,  j^(l)  =  oo.  Consequently,  d  is  not  Lipschitz 
continuous  at  s  =  1  for  small  m. 

there  is  no  net  movement  of  the  fluid  mixture  and  buoy¬ 
ancy  and  capillary  pressure  are  the  only  driving  forces.  The 
case  with  gx  =  0  is  called  capillary  redistribution.  For  these 
physical  constants  with  gx  =  —1,  the  nonlinear  coefficients 
for  two-phase  flow  are  plotted  in  figures  1  and  2  for  large 
and  small  m.  The  advective  flux  in  figure  1  is  zero  at  se  =  0 
and  se  =  1,  but  it  has  unbounded  slope  as  se  approaches 
one.  Likewise,  the  diffusion  coefficient  in  figure  2  is  zero  at 


Table  1.  Summary  of  model  parameters  for  0  <  m  <  1 

P  Q  b0 

Two-phase  3/2  +  m  1/2  +  2m  -1 
Richards  1-m  m  1 


S 


Figure  3.  advective  flux,  /,  versus  saturation,  s,  for 
Richards’  equation.  For  all  m,  /( 1)  =  —1,  and  ^(1)  = 
—oo.  Consequently,  f  is  not  Lipschitz  continuous  at 
s  =  1. 


S 


Figure  4.  Diffusion  coefficient,  d  vs.  s  for  Richards’ 
Equation.  For  all  m  d(l)  =  oo  and  ^(1)  =  oo.  Thus,  d 
is  neither  bounded  nor  Lipschitz  continuous  with  respect 
to  s. 

se  =  0  and  se  =  1,  but,  for  small  m,  the  slope  is  unbounded 
as  Se  approaches  one.  The  coefficients  for  Richards’  equa¬ 
tion  are  plotted  in  figures  3  and  4.  The  diffusion  coefficient 
itself  is  unbounded  for  Richards’  equation.  The  advective 
flux  is  finite  but  has  unbounded  slope  for  all  values  of  m  as 
se  approaches  1. 

2.6.  Asymptotic  Approximations 

We  would  like  to  understand  how  the  unbounded  slope 
of  /  and  d,  as  well  as  the  unboundedness  of  d  in  Richards’ 
equation,  affect  the  solution  behavior,  and  the  speed  of  prop¬ 
agation,  in  particular.  Since  this  behavior  occurs  only  at 
Se  =  1,  we  will  use  asymptotic  approximations  to  d  and  / 
near  se  =  1.  For  u  =  1  —  se  (i.e. ,  non- wetting  phase  satura- 
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tion)  we  have  as  u  — >  0 


i+2  m 


(20) 


dpc  m  —  11" 
ds  mm 

^  m 

kfyj  —  1  2  i 


m 


1 


krn  —  U  2 


+2  m 


—m  ,  /  —m\ 

u  +  o(u  ) 

(16) 

+  o(u  ) 

(17) 

i  /  i +2m  \ 

+  0(lL2^  ) 

(18) 

and  Richards’  equation  with 

7  ki  m  -  1  1  _ 

a  = - u 


Neglecting  constants,  we  can  approximate  two-phase  flow 
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with  the  nonlinearities 
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Making  these  approximations  yields  the  simple  model 

Ut  =  ( up)xx  +  b0(uq)x  (23) 

where  bo  =  —  1  for  two-phase  flow  and  bo  =  1  for  Richards’ 
equation.  Note  that  we  have  used  the  Kirchoff  transforma¬ 
tion  on  the  diffusion  coefficient: 


fU 

up=  d(u)du  (24) 

Jo 


u 


Figure  5.  First  derivative  of  advective  flux,  df  /dv,  ver¬ 
sus  transformed  variable  v  for  two-phase  flow.  Note  that 
df  /dv( 0)  =  —  1  in  agreement  with  the  asymptotic  approx¬ 
imation  for  f(s). 


In  terms  of  the  parameter  m,  the  approximate  two-phase 
flow  model  has  p  =  3/2  +  m  and  5  =  1/2  +  2m  where  the 
ranges  of  p  and  q  are  1  <  p  <  5/2  and  0  <  q  <  5/2.  Thus,  / 
has  infinite  slope  at  s  =  1  for  m  <  1/4,  which  holds  both  for 
the  simplified  model  and  the  original  two-phase  model.  The 
approximate  Richards’  equation  has  p  =  1  —  m  and  q  =  m 
where  the  range  of  p  and  q  is  0  <  p,  q  <  1.  Thus  |  ^  |  blows 
up  faster  than  d  when  m  <  1/2,  and  we  shall  see  that  the 
speed  of  propagation  changes  from  infinite  to  finite  at  this 
point.  For  the  van  Genuc.hten  parameter  n,  this  corresponds 
to  the  range  1  <  n  <  2,  which  is  a  case  studied  previously  in 
the  context  of  numerical  approximations  Miller  et  al.  [1998]. 
These  relations  are  summarized  in  table  1. 

2.7.  Variable  Transformation 

The  asymptotic  approximations  above  not  only  yield  sim¬ 
plified  model  equations  but  also  a  change  of  variables  for  the 
original  model  equations  that  will  yield  once-differentiable 
coefficients  over  the  entire  range  of  saturations  0  <  se  <  1. 

First,  note  that  if  we  use  the  Kirchoff  transformation 

f‘U 

(j)(u)  =  /  d(u)du  (25) 

Jo 


u 


Figure  6.  Diffusion  coefficient,  d,  vs.  transformed  vari¬ 
able  v  for  Richards’  equation,  d  is  now  bounded  and 
continuous  but  for  m  <1/2  d(0)  =  0. 


then  cj>  is  once  differentiable  as  long  as  d  is  continuous.  Hence 
only  in  the  case  of  Richards’  equation,  where  d  is  unbounded, 
must  we  consider  the  behavior  of  d. 

For  two-phase  flow  in  the  case  m  <  1/4,  we  will  make  the 
substitution 

i  2m  , 

v  =  —  «5+2m  (26) 

m 

This  transformation  yields  nonlinear  coefficients  s,  /,  and 
( f>  that  are  differentiable  in  v.  Figure  5  shows  df  /dv.  For 
m  >  1/4  no  variable  transformation  is  necessary  so  we  set 
v  =  u. 

The  variable  transformations  are  slightly  more  complex 
for  Richards’  equation.  For  m  >  1/2,  d  — >  oo  and  d/|  — > 
oo  as  u  — >  0.  We  design  a  change  of  variables  v(u)  so  that 
(dpc/ds)(dv/du)( 0)  =  —  1.  Such  a  v(u)  will  yield  a  diffusion 
coefficient  that  is  bounded  and  continuous  and,  therefore,  a 
potentiai  (f>  that  is  differentiable  and  an  advective  flux  that 
is  zero  at  u  =  0.  This  consideration  leads  to 

m  1  —  m 

v  =  m  u  (27) 

For  m  <  1/2,  d/|^|  — >  0  as  u  — >  0.  Thus,  as  we  did  in  the 
two-phase  case,  we  construct  a  change  of  variables  so  that 
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the  advective  flux  is  differentiable: 


where 


1  m 

v  =  2-  um  (28) 

m 

For  the  simplified  model  equation  we  can  use 

v  =  (29) 

These  variable  transformations  do  not  change  the  charac¬ 
teristic  speeds  for  the  first-order  part  of  the  equation,  which 
are 


dx  =  J£  =  df_ 

dt  ds 


(30) 


The  shock  speeds  are  also  unchanged,  since  for  a  disconti¬ 
nuity  with  left  state  S-  =  s(v~)  and  right  state  s+  =  s(ii+) 
we  have 


d=f(S+)-f(S-)  =  f(u+)-f(u  ) 

S+  —  S-  s(u+)  —  s(U-) 

We  make  these  observations  since,  in  the  case  of  first-order 
conservation  laws,  not  all  transformations  of  the  equations 
yield  the  same  solution  in  the  sense  of  entropy-satisfying 
weak  solutions.  In  such  cases  two  versions  of  an  equation 
that  are  equivalent  for  classical  solutions  nevertheless  have 
shock  waves  moving  at  different  speeds  for  the  same  left 
and  right  states  (c.f.  [Smoller,  1991]).  The  simple  vari¬ 
able  transformations  above  yield  equations  with  the  same 
entropy-satisfying  weak  solutions  as  the  original  equation. 
On  the  other  hand,  the  diffusion  coefficient  (in  the  sense  of 
dfi/dv)  can  change  dramatically  in  the  transformed  equa¬ 
tion,  and  the  way  it  changes  yields  additional  insight  into 
the  character  of  solutions.  For  Richards’  equation  the  dif¬ 
fusion  coefficient  of  the  transformed  equation  (m  <  1/2)  is 
now  zero  at  s  =  1  (u  =  0),  as  is  shown  in  figure  6. 

In  the  next  section  we  apply  the  results  of  Gilding  and 
Kersner  [1996]  to  characterize  the  models  above  with  re¬ 
spect  to  the  speed  of  propagation. 


3.  Speed  of  Propagation 

3.1.  Simple  Model 

Conveniently,  Gilding  and  Kersner  [1996]  provide  a  corol¬ 
lary  to  their  main  results  which  describes  the  speed  of  prop¬ 
agation  of  the  simple  nonlinear  advection-diffusion-reaction 
equation 


ut  =  ( up)xx  +  boux  +  cqu  (32) 

over  a  wide  range  of  values  of  p,  q ,  r,  bo,  and  Cq.  The  rel¬ 
evant  portion  of  this  corollary  states  that  for  Co  =  0  and 
bo  ^  0,  equation  32  has  finite  speed  of  propagation  if  and 
only  if 

1.  p  >  1  and  q  >  1,  or 

2.  q  <  1,  &o  >  0  and  p  >  q 

For  the  simplified  two-phase  flow  model,  the  first  case 
holds  when  m>  1/4  while  the  second  case  holds  for  gx  >  0. 
Thus  the  solution  will  display  infinite  speed  of  propagation 
when  gx  =  —1  and  m  <  1/4,  in  other  words,  in  the  counter- 
current  flow  problem  when  the  advective  flux  has  infinite 
slope  at  se  =  1.  For  the  simplified  Richards’  equation  model, 
only  the  second  case  holds  and  then  only  when  m  <  1/2. 

3.2.  Two-Phase  Flow  and  Richards’  Equation 


f-U 

a(u )  =  /  d{  1  —  w)dw  (34) 

Jo 

The  formula  for  d  is  given  by  equation  5  or  equation  8.  For 
two-phase  flow  we  have  b  =  f(l—u)  where  /  is  given  by  equa¬ 
tion  6.  For  Richards’  equation  we  take  b  =  f(l  —  u)  +  l  where 
/  is  given  by  equation  7.  With  these  definitions  equation  33 
is  equivalent  to  equation  4  and  a(0)  =  6(0)  =  0,  which  is 
the  form  of  the  equation  needed  to  apply  the  following  two 
results  from  Gilding  and  Kersner  [1996]. 

Theorem  1  (Gilding  and  Kersner)  The  equation 

ut  =  [a{u)]xx  (35) 

displays  finite  speed  of  propagation  if  and  only  if 

fS  1 

/  —da(w)  <  oo  for  some  5  €  (0,oo)  (36) 

J  0 

where  the  integral  is  in  the  sense  of  Lcbesgue-Stieltjes. 
Since  a  is  differentiable  on  (0,  <5),  da(w)  =  d(  1  —  w)dw.  For 
two-phase  flow  we  see  from  equation  19  that  equation  36 
reduces  to 


Cw%+m  +  o(w5+m) 
w 


dw  <  oo 


V<5  >  0 


(37) 


so  that  in  the  case  of  capillary  redistribution  two-pliase  flow 
has  finite  speed  of  propagation  for  all  m.  For  Richards’ 
equation  we  see  from  equation  21  that  equation  36  reduces 
to 


Cw~m  +  o(w~m) 
w 


dw  =  oo 


V<5  >  0 


(38) 


so  that  in  the  case  of  capillary  redistribution,  Richards’ 
equation  has  infinite  speed  of  propagation  for  all  m. 

When  /  0  we  have 

Theorem  2  (Gilding  and  Kersner)  The  equation 


ut  =  [a(u)x  +  b{u)\x  =  0  (39) 

displays  finite  speed  of  propagation  if  and  only  if 

max  {— b(w),  0}  =  0(w)  as  w  J.  0  (40) 


and 

f5  1 

/  - z—, — r - Tda(w )  <  oo  for  some  <5  €  (0,oo)  (41) 

J0  max  {b(w),w} 

For  two-phase  flow  we  see  from  equation  20  that  equation 
40  fails  when  m  <  1/4.  When  m  >  1/4  then  equation  41 
reduces  to  equation  37  so  that  two-phase  flow  displays  finite 
speed  of  propagation  in  the  gravity  segregation  case  if  and 
only  if  m  >  1/4. 

For  Richards’  equation,  equation  40  holds  for  0  <  m  <  1 
while  equation  41  reduces  to 

Kw~2m  +  o{w~2m)dw  (42) 


Table  2.  Summary  of  Finite  Speed  of  Propagation  Results 


Since  we  are  considering  the  behavior  near  se 
consider  the  equation  in  terms  of  u  =  1  —  se 

ut  =  [a{u)x  +  b(u)]x 


^  we  Capillary  Redistribution  Gravity  Segregation 

Two-phase  0<m<l  1/4  <  m  <  1 

(33)  Richards  —  0<m<l/2 
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Hence  Richards’  equation  displays  finite  speed  of  propaga¬ 
tion  in  the  gravity  segregation  case  if  and  only  if  m  <  1/2. 

These  results  are  in  agreement  with  those  for  the  simple 
model  and  are  summarized  in  table  2. 

4.  Analytical  Solutions 

In  the  following  sections  we  examine  solutions  with  and 
without  finite  speed  of  propagation  by  looking  at  special 
solutions. 

4.1.  First-Order  Approximation 

First  we  consider  solutions  of  the  first-order  approxima¬ 
tion  of  equation  4 


St  +  fx  —  0 


(43) 


on  the  spatial  domain  (— oo,  oo)  given  Riemann  initial  data, 


s( x,  0) 


s_  for  x  <  0 
s+  for  x>d 


(44) 


Equation  43  is  usually  referred  to  as  the  Buckley-Leverett 
equation  when  /  is  given  by  equation  6.  This  equation  gen¬ 
erally  only  has  weak  solutions,  which  may  be  discontinuous 
and  which  furthermore  are  not  unique  (c.f.  [Evans,  1998; 
Leveque,  1990]).  When  the  solution  fails  to  be  continuous, 
a  discontinuity  with  left  state  s_  and  right  state  s+  must 
propagate  with  speed 


f(s+)  -  f(s-) 

S+  —  S- 


(45) 


Unique  weak  solutions  can  be  specified  by  requiring  that 
discontinuous  solutions  be  in  some  sense  a  limit  of  solutions 
of  equation  43  with  a  second-order  regularization  term  euxx 
(a  vanishing  viscosity  solution).  For  non-convex  /,  this  re¬ 
quirement  can  be  stated  as  the  Oleinik  entropy  condition 
[Oleinik,  1957]  (c.f.  [Osher,  1984;  Leveque,  1990]) 

/(«)-/(«-)  /(*)-/(«+)  (46) 

S  S—  s  -  s+ 


Figure  7.  Riemann  solution  for  the  Buckley-Leverett 
equation  for  the  cases  m  =  5/6  and  m  =  1/6. 


For  example,  if  we  consider  values  of  s_  close  to  s+  we  can 
produce  a  shock  wave  moving  arbitrarily  fast.  Furthermore, 
if  we  consider  the  initial  data  reversed  (also  known  as  capil¬ 
lary  rise)  we  again  obtain  a  rarefaction  wave  with  values  of 
=  1  moving  left  at  infinite  speed:  the  saturated  region 
desaturates  instantaneously. 

For  completeness  we  note  that  the  solution  in  figure  ?? 
was  computed  automatically  by  solving  the  optimization 
problem  [Osher,  1984] 

s(0  =  argmins6[3_  s+]  [f(s)  -  £s]  (48) 

where  £  =  x/t.  This  formula  is  a  concise  statement  of  the 
entropy  satisfying  weak  solution  to  the  Riemann  problem 
when  s+  >  s_  (a  similar  formula  holds  for  s+  <  s_).  This 
problem  was  solved  in  MATLAB  using  the  fminbnd  routine 
with  tolerances  of  l.Oe  —  5. 


for  all  s  between  s_  and  s+.  For  a  general  Riemann  problem 
with  S-  <  s+,  the  unique  solution  can  be  determined  from 
the  convex  hull  of  the  set 

{(s,y)\s-  <  s  <  s+  and  y  >  f(s)}  (47) 

The  convex  hull  will  consist  of  the  graph  of  /  and  some  set  of 
chords  enclosing  regions  where  /  is  concave  down.  A  chord 
on  the  boundary  of  the  convex  hull  represents  two  values 
of  the  solution  connected  by  a  shock  wave,  and  the  slope 
of  the  chord  is  the  shock  speed  c.  We  are  interested  in  the 
counter-current  flow  problem  where  s+  =  1.  For  two-phase 
flow  we  see  from  figure  1  that  the  solution  will  have  at  least 
one  shock  wave  connecting  to  s+  =  1  and  moving  to  the 
right  (c  >  0)  when  m  is  large,  specifically  when  m  >  1/4. 
On  the  other  hand,  when  m  is  small,  the  solution  near  s+ 
will  be  a  rarefaction  wave  with  the  value  s+  =  1  traveling 
at  infinite  speed  to  the  right.  Thus,  in  this  case,  the  loss 
of  finite  speed  of  propagation  occurs  even  in  the  first-order 
approximation. 

Solutions  for  s+  =  1  and  s_  =  0  are  given  in  figure  ?? 
corresponding  to  each  case.  As  can  be  seen  from  figure  3, 
Richards  equation  always  has  solutions  consisting  of  a  single 
shock  traveling  to  the  left  because  the  flux  is  convex  down 
for  all  m.  The  solution  for  s+  =  1  and  s_  =  0  is  simply  a 
shock  wave  with  c  =  —  1.  Although  we  only  consider  this 
case  at  length,  we  can  briefly  mention  a  few  other  cases. 


4.2.  Traveling  Wave  Solutions 

There  are  two  reasons  why  studying  the  first-order  prob¬ 
lem  above  is  not  sufficient  for  understanding  the  solution 
dynamics  of  the  full  second-order  model.  The  first  is  that  a 
second-order  term  can  destroy  the  finite  speed  of  propaga¬ 
tion  property.  The  second  is  that,  even  without  a  first-order 
term,  the  equation  can  exhibit  finite  speed  of  propagation. 
In  this  section  we  will  see  how  the  full  second-order  model 
has  finite  speed  of  propagation  (in  some  cases)  by  examining 
traveling  wave  solutions. 

The  solutions  in  the  previous  section  are  the  unique  en¬ 
tropy  satisfying  weak  solutions  composed  of  continuous  rar¬ 
efaction  waves  and  discontinuous  shock  waves.  The  former 
condition  arises  due  to  the  conservation  form  of  the  equa¬ 
tion  while  the  latter  can  be  derived,  as  mentioned  above,  by 
requiring  that  any  moving  discontinuity  correspond  to  trav¬ 
eling  wave  solutions  of  a  second-order  regularization  of  the 
equation.  In  short,  we  should  already  expect  that  the  only 
traveling  waves  that  arise  in  the  full  model  are  the  shock 
waves  that  arise  in  the  simplified  model.  In  order  for  the  so¬ 
lutions  to  exhibit  finite  speed  of  propagation,  the  front  £(t) 
must  be  finite  for  0  <  t  <  r. 

To  derive  traveling  wave  solutions,  we  assume  that  the 
solution  has  the  form  v(x,t)  =  v{x  —  ct )  =  v(£)  where  c 
is  the  wave  speed  (to  be  determined).  Using  the  variable 
transformations  above,  we  can  consider  s,f,  and  <f>  to  be  C1 
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functions  of  v.  Making  this  substitution  into  Equation  4 
yields 


-CS£  +  /e  -  foe  =  0  (49) 

Integrating  once  with  respect  to  £  and  rearranging  yields 

fo  =  -ca(0  +  /(C)  +  6  (50) 

where  b  is  the  constant  of  integration.  We  will  enforce  the 
following  asymptotic  boundary  conditions  on  equation  50 


lim  v  =  v±  (51) 

£ — >±oo 

lim  fo  =  0  (52) 

£ — *■  i  OO 

and  write  f±  for  f(v±)  and  s±  for  s(v±).  These  boundary 
conditions  require 


lim  —  cs  +  /  +  b  =  0 

V — ►  U-J- 


which  implies 


b  —  cs-  -  f-  =  cs+  -  /+ 
From  equation  54  we  can  also  determine  c 

- 

s+  —  S- 


(53) 


(54) 


(55) 


Note  that  equation  55  is  the  Rankine-Hugoniot  condition. 
We  will  take  b  =  cs+  —  /+.  The  traveling  wave  equation  is 
then 

fotfo  =  f(v)  ~  |  ls(v)  ~  a+]  +  /+ j  (56) 

=  h(v-,v+,v)  (57) 


Equation  56,  with  boundary  conditions  given  by  equations 
51  and  52,  may  not  have  a  solution,  and,  if  it  does,  solutions 
are  not  unique  (if  u(C)  is  a  solution,  then  v(C(£  — fo))  is  also 
a  solution).  For  two-phase  flow  and  Richards’  equation  and 
the  variable  transforms  we  have  defined,  fou)  is  a  monotone 
decreasing  function  of  v.  Since  we  consider  the  case  v+  =  0 
and  0  <  V-  <  1,  this  implies  that  fov _)  <  (j>{v+ ).  The 
asymptotic  boundary  conditions  can  only  be  met  if  h  >  0 
for  v+  <  v  <  V-.  The  term  in  braces  on  the  right-hand 
side  of  equation  56  is  an  expression  for  the  chord  connecting 
/+  and  /_.  The  requirement  that  h  >  0  can  then  be  inter¬ 
preted  as  a  constraint  that  the  chord  connecting  /+  and  /_ 
lie  below  the  graph  of  f(s).  Rearranging  the  requirement 
that  h  >  0  yields  half  the  Oleinik  entropy  inequality  above. 
The  other  half  comes  from  using  the  second  definition  of  the 
constant  b  above.  Thus,  the  analysis  of  the  first-order  case 
carries  over  to  the  full  model  in  the  following  way:  traveling 
wave  solutions  exist  whenever  the  solution  of  the  Riemann 
problem  consists  of  a  single  shock  wave.  Traveling  wave  so¬ 
lutions  exist  for  the  two-phase  model  when  m  >  1/4  and  in 
all  cases  for  the  Richards’  equation  model. 

The  question  of  whether  the  support  of  v(x ,  t)  has  an  up¬ 
per  bound  for  t  >  0  (i.e.,  that  the  solution  has  finite  speed 
of  propagation)  can  then  be  answered  by  determining  when 
the  solution  of  equation  56  with  initial  condition  V-  >  0 
goes  to  v+  for  £  <  oo.  Dividing  equation  56  by  its  (positive) 
right  hand  side  and  integrating  once  yields 


where  Vo  =  v(0)  can  be  chosen  arbitrarily  in  (u+,  v~).  Thus, 
if  the  integral  in  equation  58  is  bounded  as  v  — ►  v+,  v(£)  goes 
to  in  “finite  £”  and  the  front  £( t )  is  given  by  tG(u+). 

Consider  first  the  simplified  model  problem  with  u_  =  1, 
u+  =  0.  The  integral  above  is 


e  = 


pu 


p- 1 


—bouq  —  cu 


-du 


(59) 


If  bo  =  —  1,  p  =  1  and  q  =  2,  the  simplified  model  is  viscous 
Burgers’  equation,  and  we  have 


£(«)  =  log(l  -  u)  -  log(u)  (60) 

Since  £( u )  is  unbounded  as  u  — >  u±,  we  see  that  the  solu¬ 
tion  has  infinite  speed  of  propagation.  On  the  other  hand, 
for  p  —  2,  we  have 


£(w)  =  21og(l-u)  (61) 


Figure  8.  Traveling  wave  solution  of  the  simplified  two- 
phase  flow  model. 


G(v) 


d<p 

dv 


f(v)  -  [c(s(u)  -  s+)  -  /+] 


dv  =  £  (58) 


Figure  9.  Traveling  wave  solutions  of  the  simplified 
Richards’  equation  model.  Note  that  the  wave  speed  is 
negative,  and  for  m  >  1/2  the  support  is  unbounded. 
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and  therefore  £(f)  =  2 1.  Note  also  that  ^  =  —  2/(1  —  u) 
and  therefore  —1/2  <  <  0:  solutions  are  Lipschitz  con¬ 
tinuous.  For  p  =  3  we  have 

£(m)  =  3(u  +  log(l  —  u))  (62) 


so  that  £(t)  =  3 1  and  the  speed  of  propagation  is  again  finite. 
However  ^  =  3[1  —  1/(1  —  u)]  and  therefore  — >  — oo  as 
x  — >  C(f):  solutions  fail  to  be  Lipschitz  continuous  in  space 
or  time  at  </(t). 

For  two-phase  flow  we  know  q  >  1  if  m  >  1/4,  so  we  can 
write  the  inverse  of  the  solution  as 

/U 

- — ,  P  . - -ydu  (63) 

U2-P(U9-!  -  1)  V  ' 

Since  p  >  1  this  implies  that  the  front  for  the  traveling  wave 
solutions  to  two-phase  flow  move  with  finite  speed  for  all 
cases  with  m  >  1/4,  in  agreement  with  the  theory. 

For  the  full  two-phase  model  there  is  no  traveling  wave 
solution  for  the  case  when  m  <  1/4.  For  m  >  1/4  the  travel¬ 
ing  wave  solution  is  bounded,  and  for  m  >  1/2  the  derivative 
of  the  solution  blows  up  at  the  front.  We  can  consider  equa¬ 
tions  61  and  62  as  approximate  models  for  the  behavior  of 
two-phase  flow  near  s  =  1.  The  traveling  wave  solution  for 
m  =  3/4  is  given  in  figure  8  (£  is  shifted  and  scaled). 

For  Richards’  equation,  on  the  other  hand,  we  have 
0  <  p,q  <  1,  p  =  1  —  q  (q  —  m),  and  we  can  see  by  re¬ 
arranging  the  integrand  of  equation  59  as 


£  = 


_ — _ du 

m2(1-p)(1  _  ul-g)aU 
1  —  m 

u2m(  1  -  u1-"1) 


(64) 

(65) 


that  we  have  infinite  speed  of  propagation  when  m  >  1/2. 
When  m  =  1/2  we  have 

£(«)  =  log(l  -  Vu)  -  logty^u));  (66) 

When  m  =  1/4  we  have 

£(«)  =  2 uV2  +  ^4  “  l0S(wl/4)  +  ^  ~  (67) 


When  m  =  3/4  we  have 


£(u)  =  log(l  —  u1^4)  —  log^1/2  +  u1^4  +  l)/2 

-t-31^2  arctan(l/3(2u1^4  +  ^S1^2)  (68) 


Traveling  wave  solutions  for  m  =  1/4  and  m  —  3/4  are  given 
in  figure  9  (£  is  shifted  and  scaled). 


4.3.  Barenblatt’s  Solution 

The  traveling  waves  in  the  previous  section  displayed  both 
finite  and  infinite  speed  of  propagation.  As  was  seen  in  the 
calculations,  the  boundedness  of  the  traveling  wave  depends 
mainly  on  the  ratio  of  diffusion  and  advection.  When  ad- 
vection  is  absent,  the  equation  can  still  have  finite  speed  of 
propagation  and  therefore  moving  fronts.  We  now  consider 
the  simplified  model  given  in  equation  23  with  g  =  0  and 
qt  =  0  (i.e.  no  advection): 

ut  -  (up)xx  =  0  (69) 

By  seeking  a  solution  of  the  form  u(x,t)  =  t~av{xt~d)  it 
can  be  shown  that  a  solution  on  (—00,00)  is  given  by  (c.f. 


[Evans,  1998]) 


u(x,  t) 


1 


tp+ 1 


b  — 


P~  1 
2  p(p  +  1) 


1*1  <  C (t) 
(70) 


for  p  /  1,  where  b  is  determined  from  the  particular  initial 
conditions  u(a:,0)  =  M5(x),  and  where  5  is  the  Dirac  delta 
function.  For  p  >  1  the  front  location  is  given  by 


Mp  + 1)  A 
p- 1  ) 


1/2 


tp+ 1 


(71) 


Thus  the  propagation  speed  is  finite  for  p  >  1,  for  all  t.  Note 
also  that  |ut|  and  \ux\  blow  up  at  £(t )  if  p  >  2.  The  case 
p  >  2  corresponds  to  two-phase  flow  with  m  >  1/2.  The 
case  where  p  <  1  corresponds  to  Richards’  equation.  Since 
in  this  case  p  —  1  <  0,  we  have  that  1/  =  00  and  the  speed  of 
propagation  is  infinite  (the  solution  is  positive  for  all  x  for 
t  >  0).  Graphs  of  the  solution  are  given  in  figure  10. 


5.  Numerical  Solutions 

We  are  capable  of  writing  closed-form  solutions  in  only 
a  few  of  the  instances  above.  We  were  able  to  demonstrate 
the  effect  of  infinite  speed  of  propagation  in  a  few  special 
cases  of  approximate  model  equations.  To  look  at  the  more 
general  cases  we  need  to  use  numerical  solutions. 

We  will  look  at  only  simple  discretization  methods  here: 
forward  and  backward  Euler  in  time  for  fixed  time  step  and 
central  and  upwind  finite  differences  in  space.  For  the  do¬ 
main  [0, 1]  x  [0,  T]  the  discrete  solution  is  defined  over  the 
spatial  grid  Xi  =  ih,  i  =  0, . . . ,  n,  h  =  1/n  at  discrete  times 
tJ  =  jk,j  =  0, . . . ,  m,  k  =  Tim.  We  consider  first  Dirich- 
let  boundary  conditions  s(0)  =  s(l)  =  S  so  that  a  gen¬ 
eral  fully  discrete  form  of  the  equation  is,  i  =  1, . . . ,  n  —  1, 

S  —  (^o , •  ■  -  ,  Sn ) 


7  7  —  1  , 

^  ^  +~h 


<Ks*+i)  -  <P(si) 


—  Fj 


i+1/2 


(s’) 


</>«)  -  0(s*-l) 


—  Fi- l/2(s*) 


(72) 


Figure  10.  Barenblatt’s  Solution 
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First  we  choose  *  =  n  and  the  numerical  flux  F  as 


^i+i/2(s) 


f{&i)  +  f(si+ 1) 
2 


(73) 


Consider  initial  conditions  s(x,  t )  =  S+w(x,  t )  where  w(x,  t) 
is  a  small  perturbation.  Linearizations  about  S  are  given  by 

0(a)  «  <KS)  +  ^(S)(s  -  S)  =  </>(S)  +  D(s  -  S)  (74) 
f(s)  «  f(S)  +  F(s-S)  (75) 

where  ^(S1)  =  D  and  ^(S)  =  F.  Substituting  these 
linearizations  for  the  nonlinear  coefficients  yields  the  stan¬ 
dard  fully  discrete  linear  advection-diffusion  equation  for 
i  =  1, . . . ,  n  —  1 


„J  _  ( Dk  _  Fk\  J- 1  i  (-i  MU-1  I  ( Dk  ,  Ffc  \  j-1 

si  -  (T?  ~  2h  lsi+l  +  “  FF-)si  +  ITT  +  2h  >si-l 


(76) 


For  periodic  boundary  conditions  the  standard  Von  Neu¬ 
mann  stability  analysis  shows  that  the  method  is  stable  for 
time  steps  k  >  0  satisfying  (c.f.  [Strikwerda,  1989;  Hunds- 
dorfer  and  Verwer,  2003]) 


k<^=k0  (77) 


If  h  is  then  chosen  such  that 


,  .  2d 

h-wr  0 

(78) 

then 

max  |  sj 

X 

< 

(79) 

(80) 

(81) 

< 

(82) 

+<i  -  ™)Wr'\ 

(83) 

,Dk  Fk ..  ,_i,l 

+W+2ft)^-'} 

(84) 

< 

(85) 

X 


and,  therefore,  saturations  will  remain  in  the  physical  range 

0  <  s  <  1. 

From  equations  19-22  we  have  for  two-phase  flow 


k0 

h0 


2m  2' 


(|  +  2  m)u 


2m—  ^ 


and  for  Richards’  equation 


ko 

ho 


h2 

2u~m 

u~ 


(1  —  m)urn~1 


(86) 

(87) 


(88) 

(89) 


which  yields  the  following  three  cases: 

1.  (Two-phase  Flow)  ho  — >  0  as  S  — >  1. 


2.  (Richards’  equation  for  m  >  1/2)  ko  — >  0  as  S  — >  1 
but  ho  is  bounded  away  from  zero. 

3.  (Richards’  equation  with  m  <  1/2)  ho  — >  0  as  S  — >  1. 
Since  fco  — >  0  as  ho  — >  0  the  discretization  is  essentially 

unusable  for  these  model  equations.  If  we  use  an  upwind 
advection  discretization  given  by 


Fi+  l/2(s) 


Fsi  if  F  >  0 
Fsi+ 1  if  F  <  0 


(90) 


then  the  discrete  maximum  principle  holds  [Strikwerda, 
1989]  if 


2 Dk  Fk  „  „ 


(91) 


This  condition  can  be  met  only  in  the  two-phase  flow  case 
with  m  >  1/4.  We  therefore  elect  to  use  backward  Euler 
with  upwinding: 


J  —  J-1  I  DkJ  __  (2Dk  _l  Ek\J  _i_  (UF  _i_  Ffc \  j 

Si  —  si  T  h2  si+1  I  h 2  T-  ft.  JSi  +  l  T  h  )Si_  1 


(92) 


Figure  11.  Solutions  of  the  simplified  two-phase  model 
at  t  =  0.05. 


1 1— 
0.8  - 
0.6  - 

0.4  - 

0.2  - 

0  — 
0 


- p=3,q=3 

-  -  p=3/2,q=1/2 


Figure  12.  Solutions  of  the  simplified  two-phase  model 
at  t  =  0.5. 
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Taking  absolute  values  and  assuming  F  >  0  yields 


2  Dk  Fk . .  j. 

<I  +  ^  +  T  )«l 


Dk  ,  ,  Dk  Fk .  , 

+  T2-*i+i  +  (-p- +  —  K-il 


< 


h 2 

kn+fv 


i+ll 


.  Dk  Fk . ,  a 

+(^  +  x)|s^ 


(93) 

(94) 


Thus  we  have  unconditionally 


max  I  si  I  <  max  I  si  1 1 


(95) 


Since  backward  Euler  is  fully  implicit,  we  must  solve  a  sys¬ 
tem  of  nonlinear  algebraic  equations  at  each  time  step  of 
the  form  G(s3)  =  0.  Standard  theories  for  the  convergence 
of  Newton’s  method  do  not  apply  when  G  is  not  Lipschitz 
continuous,  which  is  the  case  when  m  <  1/4  for  two-phase 
flow  or  m  <  1/2  for  Richards’  equation  [ Kelley ,  1995;  Qi 
and  Sun,  1993].  We  therefore  apply  the  change  of  variables 
as  described  above  to  solve  the  equivalent  root-finding  prob¬ 
lem  G(vJ).  This  equivalent  problem  is  Lipschitz  continuous 
in  v,  and,  therefore,  the  convergence  theory  for  non-smooth 
Newton  methods  is  applicable. 

In  computations  with  the  two-phase  and  Richards’  mod¬ 
els,  <j>(v)  is  approximated  with  a  piecewise  linear  spline.  To 
compute  the  spline  approximation  we  set  vi  —  lAv  where 
Av  —  ma xv /(ns  —  1)  and  ns  is  the  number  of  spline  knots. 
We  then  set  (po  =  0  and  approximate  <p  recursively  as 


(pi  =  tpi-!  +  A vd(vi)  l  =  1, . . . ,  ns  (96) 


In  the  nonlinear  context  we  use  the  following  simple  gener¬ 
alization  of  the  upwind  formula 


^i+i/2(s) 


f(Si) 
/(««+ 1) 


if 

if 


fi  +  l—fi 

si  +  l~si 

fi+l—fi 

Si  +  l—Si 


>  0 
<  0 


(97) 


This  flux  is  not  generally  appropriate  for  the  first-order  ap¬ 
proximation  to  two-phase  flow  (Buckley-Leverett),  since  the 
flux  is  non-convex.  It  is  not  entropy  satisfying  in  the  case  of 
transonic  rarefaction  waves  [Leveque,  1990]. 

The  resulting  nonlinear  system  was  solved  with  Newton’s 
method  to  a  relative  residual  tolerance  of  10-5.  Figures  11 
14  are  plots  of  the  solution  to  the  simple  model  for  param¬ 
eters  corresponding  to  two-phase  flow  and  Richards’  equa¬ 
tion.  Figures  15-16  are  plots  of  the  solution  to  the  two-phase 
flow  model  at  t  =  0.05  and  t  =  0.5  m  =  1/5  and  m  =  5/6. 
The  consequence  of  infinite  speed  of  propagation  is  the  thin 
wedge  of  saturation  that  extends  all  the  way  to  the  right 
end  of  the  domain  as  opposed  to  the  sharp  moving  bound¬ 
ary  that  forms  in  the  finite-speed  case.  Similar  behavior  is 
seen  in  the  solutions  of  Richards’  equation  shown  in  figures 
17-18.  The  sharp  boundary  traveling  to  the  left  in  both 
two-phase  flow  and  Richards’  equation  is  a  consequence  of 
finite  speed  of  propagation  of  wetting  phase  into  non- wetting 
phase  (the  degeneracy  at  se  =  0). 


6.  Conclusions 

We  have  considered  two  models  of  two-phase  flow  and  the 
the  closure  relations  presented  in  Parker  et  al.  [1987],  which 
are  based  on  the  relative  permeability  model  of  Mualem 
1976]  and  the  capillary  pressure  model  of  van  Genuchten 
1980].  For  horizontal  capillary  redistribution  (zero  grav¬ 
ity  and  zero  total  fluid  velocity),  the  full  two-phase  model 
displays  finite  speed  of  propagation.  Richards’  equation  dis¬ 
plays  infinite  speed  of  propagation  in  this  case. 

For  two-phase  flow  in  the  case  of  counter-current  flow 
with  the  pore  size  distribution  index  m  <  1/4  (n  <  4/3), 
non-wetting  phase  propagates  into  wetting  phase-saturated 


regions  with  infinite  speed.  In  this  case,  infinite  speed  of 
propagation  is  a  consequence  of  the  unbounded  derivative 
of  the  advective  flux.  At  sw  =  1,  the  diffusion  coefficient  is 
degenerate  while  the  characteristic  speeds  of  the  underlying 
first-order  equation  are  infinite.  A  standard  upwind  advec- 
tion  discretization  with  forward  Euler  time  stepping  fails  to 
have  a  discrete  maximum  principle  in  this  case,  while  back¬ 
ward  Euler  with  upwinding  has  no  such  requirement. 

For  Richards’  equation,  in  the  case  of  counter-current  flow 
with  m  <  1/2  (n  <  4/3),  air  propagates  with  finite  speed. 
As  in  the  degenerate  case  of  two-phase  flow,  however,  the 
standard  upwind  advection  discretization  with  forward  Eu¬ 
ler  time  stepping  fails  to  have  a  discrete  maximum  principle. 
Again,  the  backward  Euler  method  has  no  such  requirement. 

Non-Lipschitz  continuous  closure  relations  for  two-phase 
flow  and  Richards’  equation  both  cause  non-trivial  differ¬ 
ences  in  the  qualitative  behavior  of  solutions  and  in  the 
difficulty  of  obtaining  numerical  solutions.  While  variable 
transformations  can  be  used  to  derive  equivalent  models 
with  smooth  coefficients,  the  most  important  numerical  dif¬ 
ficulties  arise  from  the  properties  of  the  solution  itself  and 
are  therefore  unaffected  by  the  reformulation. 
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Figure  13.  Solutions  of  the  simplified  Richards’  model 
at  t  —  0.05. 


X 


Figure  14.  Solutions  of  the  simplified  Richards’  model 
at  t  =  0.5. 


Figure  15.  Solutions  of  the  two-phase  model  for  m  = 
5/6  and  m  =  1/6  at  t  =  0.05.  The  loss  of  finite  speed 
of  propagation  for  the  m  =  1/6  case  causes  the  entire 
domain  to  become  unsaturated  instantly. 


Figure  16.  Solutions  of  the  two-phase  model  for  m  = 
5/6  and  m  =  1/6  at  t  =  0.5. 
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Figure  17.  Solutions  of  Richards’  equation  for  m  =  5/6 
and  m  =  1/6  at  f  =  0.05.  The  loss  of  finite  speed  of  prop¬ 
agation  for  the  m  =  5/6  case  causes  the  entire  domain  to 
become  unsaturated  instantly. 


Figure  18.  Solutions  of  Richards’  equation  for  m  =  5/6 
and  m  =  l/6att  =  0.5. 


